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A 

AA 

{a} 

a.0 “3.7 
B 

BB 

C 

Cl 

c 

c2 

D 

DD 

DD®P 

E 

E 

[F] 

{f} 

G 
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[K] 

[K]^P 

L 

Lw 

Le 

[M] 

{M} 


Part of elemental stiffness matrix 
Slope of stress strain curve 
Flow vector 

Constant used in Newmark's Implicit scheme 
Part of elemental stiffness matrix 
Matrix relating strains and displacements 
A part of elemental stiffness matrix 
Property matrix 
Dilatational wave velocity 
Shear wave velocity 
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Stress strain matrix 
Elastoplastic stress strain matrix 
A part of stress strain matrix 
Modulous of elasticity 
A part of elemental mass matrix 
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Elastic thickness 
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Elastoplastic stiffness matrix 
Differential operator 
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A part of elemental stiffness matrix 

Uniformly distributed load 

Shear force 

Edge moment 
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Displacement in X-direction 

Displacement in Y-direction 

Displacement in Z-direction 

Constants used in Newmark's Implicit scheme 

Strain 

Boundary- 

Shear strain 

Poisson's ratio 

Stress 

Rotation 

Area 
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OSTOCTT 


In practical cases whet a pLate is subjected to a dynamic 
loading , the geometric: and ia-ater~ia.l linearity may not be 
preserved. This study is aimed to investigate the effect of 
plasticity in the dynamic response . _As there are no closed form 
solutions available owing to the complexity of the problem 
, numerical methods are used to s oL -ves the problem. 

For this computational work an efficient FEM code is 
developed using the "robust trdang~ul ar element" developed by 
Zienkiewiz [1991] . "Mindlin pLate ttieory" in the form of mixed 
formulation is used for modelling pLate bending. The inertia 
effect of the element is repre s&ntedlloy consistent mass matrix. 

As the stresses are var^^ing along the thickness, the 
plasticity should be applied together with elasticity for every 
iteration. This introduces the concept of "elastic height" , which 
decides the stresses in the plate. And to find out the 
elastoplastic s tiff ness matr ic es., th_e integration should be done 
in both elastic core and plastic couter region. 

The computer code is first -used, to solve a standard problem 
and after validation it was ns ed Eoi: dynamic cases of impact 
loading. The dynamic loading iincxeas«es the plasticity and in turn 
the plasticity in the material acts ss an dampener to the dynamic 
oscillatory response . The iaterdlepeoclence of plastic and dynamic 
behaviours is analyzed in ttie present work. 



1.1 GENERAL INTRODUCTION 


Plasticity in the plate is different from a simple two 
dimension analysis of the plasticity. The main difference is the 
variation of stress along the thickness. As the first 
approximation of the plate stress and strain in the thickness 
direction are neglected [Reissner Mindlin plate formulation] . And 
if the plate is thin then the out of plane stresses are 

also neglected. For a thick plate in-plane-stresses 
varies with the height while out of plane stresses remain almost 
constant along the height. This variation of the stresses affects 
the size of the elastic region in the core of the plate. 

In general plate formulation is made in terms of lateral 
forces and moments. In plastic case since the in-plane stresses 
are different in the elastic core and outer plastic zone, the 
calculation of the stress resultant moments becomes a function 
of the size of elastic core. The whole thickness is not plastic, 
proper weightages are to be given to the elastic and plastic 
stiffness matrices to calculate effective stiffness matrix. 
Calculations of unbalances after the integration of dynamic 
analysis and the plastic analysis is to be done to incorporate 
both the aspects properly. These taken together increases the 
complexity of the problem. 

1.2 LITERATURE SURVEY 

The thin plate theory developed by Kirchhoff is modified for 
thick plates by Reissner and Mindlin by relaxing the assumption 
that the normals to the midplane remain normal. This was done 
since the shear force terms can not be neglected in many loading 
and support conditions. Reissner and Mindlin theories can be 
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extended to thick plates and thus known as Reissner-Mindlin 
formulation. Regarding Finite Element Method a difference 
between thin plate and thick plate analysis is that in former it 
is possible to represent the state of deformation by one quantity 
deflection w and thus c’ continuity of w is needed, while in 
latter w can have C° continuity along with rotation and shear 
degree of freedom. 

While analytical solutions remain difficult for thick 
plates, the advent of Finite Element Method made thick plate 
theory simpler to implement. Thick plate theory on the imposition 
of additional assumptions behaves as thin plate theory. It is 
more natural to start with thick plate formulation and make 
transition to thin plate theory by imposition of certain relation 
between degrees of freedom [Zienkiewicz and Taylor, 1991] . In 
thick plate formulation instead of one approximate function for 
w ( for calculation of rotational angle ) , independent 
interpolation functions for deflection w and rotation 6 are 
taken. During formulation, if the shear forces are eliminated, 
the formulation is known as "irreducible formulation" and when 
functions for shear forces are also retained the formulation is 
known as "hybrid formulation" . The latter formulation gives more 
degree of freedom and so it is preferred to work with in most of 
the cases [Zienkiewicz and Taylor, 1991] 

Another aspect is the selection of the element to perform 
satisfactorily for thick plates as well as thin plates. Many 
elements are attempted to avoid the problems of locking and 
singularity but the results of "triangular robust element" are 
found favourable in this regard [Zienkiewicz and Lefebvre, 1988] . 

For determination of plasticity three things are required: 
an yield criteria, flow rule and plastic stress strain law. A 
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fundamental difference between plastic and elastic analysis is 
that in elastic solution total stress can be evaluated from the 
total strain alone, while in the plastic response evaluation of 
stress strain at any time depends on the stress strain history. 
This depends on the experimental results which are termed as 
yield criteria. Commonly used yield criterion are continuous one 
but piecewise linear yield criterion can also be used [Giulio 
Maier et al.,1987]. 

Yield criterion defines the onset of plasticity and flow 
rule governs the relationship between stress and plastic strain 
increment after the yield point. Flow rule requires the 
definition of a potential which is taken to be the same as yield 
surface for the case of associated plasticity. The plastic stress 
strain law has been defined for nonlinear, linear and perfectly 
plastic cases [Yamada et al., 1968]. In latter two cases the 
total stress strain curve is bilinear and in former case usually 
some power law satisfying the experimental results is defined. 
Flow rule, thus can be used in plastic cases by forming a plastic 
potential and then defining the plastic strains by the same. When 
the plastic potential surface is the same as the yield surface, 
the plasticity is termed as "associated" otherwise the plasticity 
remains "unassociated" . The flow rules for loading and unloading 
are not same, and that makes the problem difficult. 

Out of many iterative schemes to solve the above non-linear 
problem, Newton-Raphson method is found favourable in terms of 
ease in mathematical formulation and in computational effort 
[Simo and Taylor, 1986] . In Newton-Raphson method, incremental 
strain is multiplied to plastic stress strain matrix to get the 
stresses, while in Runge-Kutta method the incremental strain is 
subdivided in many subincrements then individual subincremented 
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strains are multiplied to their corresponding stress strain 
matrices to get subincremented stress . These stresses are added 
up to get incremental stress. This is known as 
"subincrementation" [Nayak, 1971 ; Bathe, 1990]. 

For plasticity in the Mindlin plates, it is suggested to 
use full three dimensional analysis of stress field. Inclusion 
of shear terms in the yield function allows the spread of 
plasticity from the extreme fibre to entire plate thickness 
[Hinton and Owen, 1980] . A very thick plate can be divided in many 
layers of thin plates and with the Mindlin plate theory, the 
problem is solved [Hinton and Owen, 1980]. But the triangular 
robust element is proved favourable for transition from thick to 
thin plates, and the non layered approach can be used for 
moderately thick plates [Zienkiewicz and Lefebvre, 1988] . 
Elastoplastic numerical experiments has been carried out adopting 
a full three dimensional version and different criteria of yield 
in the past [Hinton and Owen, 1984] .These results show excellent 
agreement in static elastic analysis and the prediction of 
vibrational frequencies and buckling loads. These results forms 
the basis of comparison for the present work. 

In elastodynamic, a detailed analyses of various time 
integration schemes from point of view of accuracy and speed, 
were reported by [Seron et al.,1990; Zienkiewicz and Taylor, 1991 ; 
Bathe, 1990] . It is concluded that, even though central difference 
scheme is the best in terms of computational speed and cost, 
Newmark's constant average scheme gives better accuracy and 
offers unconditional stability. Whereas, the time step size has 
negligible effects when central difference scheme is employed 
[Bazant and Celep,1983]. 
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In elastoplastic dynamic analysis based on Mindlin plate, 
the work started quite late. Most of the work is in the area of 
concrete structure under dynamic loading [Kujawski et al.,1984] . 

Experimental information on material behaviour is scarce 
for most structural materials [Hinton and Owen, 1980], 
Instantaneous yield stress is significantly influenced by the 
rate of straining, also the value of elasticity modulus is found 
to be dependent on the strain rate for structural material those 
varies from the materials of limited ductility to brittle elastic 
behaviours. For this purpose a better understanding of obseirved 
phenomenon and underlying microscopic behaviour is required 
[Hinton and Owen, 1980] 


In loading and unloading the yield point of the material 
changes because of strain hardening. A model for yield limit 
degradation has been developed by Bazant [ 1987 ] . Other study 
[Panzeca et al.,l987] deals with a transformation matrix that 

takes care of yield point variation according to the loading and 
history. 

Application of numerical model to dynamic flow of steel 
specimen has been done by Casadei et al.[i987] . The results 
compared favourably with the experimental result. For beams and 
rods a code validation case study is done by Marques [1987] . An 
elasto-plastic one degree of freedom numerical model was 
presented for steel slab under blast loading [Mazza et al.,1987] . 
The results were compared with static elastoplastic loading. Also 
both of the results were compared with experimental data. It was 
found out that in the dynamic loading conditions the 
displacements were intensified at the plastic zone compared with 
those of static loading. In the area of explosive welding also 
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the study was carried out [Ante, 1987] . In this case welding of 
two concentric tubes were considered while explosive was put 
inside the inner tube . The specimen experienced both mechanical 
as well as thermal shock. In this problem plasticity was induced 
by thermal load also. 


PRESENT WORK 

The present work deals with the study of the effect of 
plasticity in plates under various type of loadings both static 
and dynamic . 

A mixed finite element analysis using Reissner Mindlin 
plate theory is used to model the plate bending. To avoid locking 
and singularity a triangular robust element (T6/3 B3) [Zienkiewicz 
and Taylor, 1991] is used. 

Chapter 2 deals with the theoretical part of the 
formulation of plate bending and plasticity. Chapter 3 deals with 
the finite element implementation of the plastic analysis and its 
integration with dynamic analysis. Chapter 4 describes the 
validation of the code and the behaviour of plasticity in dynamic 
loading and effect of plasticity on the dynamic response. Chapter 
5 presents the conclusions and the suggestions for future work. 
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CHlLtx'BK 2 


BASIC EQUATIONS OF ELASTOPLASTIC DYNAMIC PLATE ANALYSIS 

2.1 INTRODUCTION 


This chapter presents the basic equations governing the 
plate under elastoplastic and dynamic conditions. As explained 
in earlier chapter, by appropriate element selection and 
formulation Mindlin plate theory can be used for thick and thin 
plates . 

The Mindlin plate theory makes the following assumptions, 

1. The strain and the stress components in Z-direction are 
negligible ( cr^ and « 0 ) (fig 2.1), where Z-direction 
represents the thickness direction. 

2. The normals to the middle plane of the plate remains straight 
during deformation. 

Applying above assumptions, static ecpailibrium equations 
for plate bending problem can be obtained [Timoshenko, 1959] in 
terms of midplane rotations (6jj,0y) , lateral displacements (w) of 
the midplane and the shear f orces (S^^, S^) at the cross-section. 
Sign convention and the variables are shown in fig(2.1). 

Stress resultant moments and shear forces are defined in 
terms of stress variables as follows: 
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2.1 
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(2.1a) 


f*h/2 

z dz 

. -h/2 ^ 


'‘A/2 

Af = Oy z dz 

- A/2 


''A/2 

^xv = '^xy ^ 

^ - -h/2 ^ 


(2.1b) 


-A/2 -A/2 

^x = "^xz ' -^v = "^yz 

--A/2 ^ --A/2 


(2.1c) 


For the dynamic analysis inertia forces due to mass and 
moment of inertia are considered in equilibrium equations. 
Formulating equilibrium equation using d' Alembert's principle, 
the following equations are obtained. 


dM^ 

dx 




^^xy 

dy 



= 0 


(2.2a) 


dx 


+ 


dMy 

By 


P— 6. = 0 
^ 12 ^ 


(2.2b) 


+ q - phw =0 (2.2c) 

dx dy 


Strains at a height of z from the centre plane are given as. 
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(2.3a) 


ae. 


30X 

e. = z , e, = z 


dtJx 

ay ^ ax . 


y = 6 ^ Y.e+-^ 

^ ax ' ay 


(2.3b) 

(2.3c) 


2.2 THE MIXED FINITE ELEMENT FOKMDLATION 

The finite element formulation involves three sets of 
variables 0,S,w. These variables at any point within the element 
are expressed in terms of nodal parameters (6 0y, S^, Sy,w) of the 
element using shape functions. 


0 = 



Ngd 


(2.4a) 


5 = 



= 


(2.4b) 


w = N^,w (2.4c) 

Where 0,S,w are the nodal values of respective variables. 

Now using eg (2. 3) and eg (2.1) the eg (2. 2) can be rewritten 
as follows. 
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(2.5a) 


L ^DLe - S 



= 0 


0 - C^S + w = 0 


(2.5b) 


v^S + g - phw = 0 


(2.5c) 


Where L,V operators and C, matrices are given by. 


_d_ 

dx 


0 


L = 0 


- 


_a_ 

dy 

d 

dx 


'A.' 

dx _ 1 1 0 

d ' ■ "^0 1 

L dy. 


D 


1 V 

Eh^ V 1 


12(l-v2) 0 0 


0 

0 

1-v 

2 


Now by Galerkin method, using eq(2.5) and eq(2.4), and 
then integrating over the area of the element, the following 
equations are obtained. 


AQ + BS + Ed = fQ (2.6a) 


- PS + Cw = 0 (2. 6b) 
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(2.6c) 


C^S + Fw = f„ 


Where, 


A= ^ {LNQ)^D(LNQ)dQ , B= ' Ne^N^dQ 


C= ^N^ivNjdQ , P= ^NsC.N^dd 

vQ -Q 


fe=''Ne^Tdr , f„='Njg dCi 

•Tr -Q 


E = 


p 

12 


" Nq^ n. dCi 
o “ 


F = p h' nJ dCl 

.Q 


In the above equations the vector T represents two 
moment components imposed on the boundary by the traction and q 
is the intensity of the imposed lateral force. The shape 
functions Nq,N^^ are chosen to have Cq continuity, but because of 
shear forces are defined only for inside nodes of the element, 
N can be discontinuous between the elements . Such a choice 

S 

allows S to be defined locally for each element and thus be 
eliminated when P is a non- singular matrix. Thus the problem 
includes only 6 and w as variables . 

The final formulation in terms of 6 and w, after 
condensing S, becomes : 
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(2.7) 


A + BP-^'C 0 ^ O =.0 A 

_ + ^ _ = 

c'^p'^B^ c'^p~^cy. ^.y. \_^v 

In short eg (2. 7) can be written as, 


[ir]®{u} + [M]®{u} = {F}® 


( 2 . 8 ) 


2.3 PLASTIC ANALYSIS 


When the stresses go beyond the yield value, then stiffness 
matrix changes depending on the individual stress component, 
effective stress and flow rule. In such cases no direct solution 
can be used, so iterative schemes are used to solve the problem. 

To do the non-linear problems, the following steps are 

used: 

a. Calculation of stress and strain 

After obtaining the maximum possible elastic solution, the 
correspondig strains are calculated by the following formula 
which is obtained by substituting eg(2.3) using the variables 
from eg (2 .4) . 


{e} 


z{LNq) 

Ne 


0 

{u} = {BB] {u} 

vN„_ 


(2.9) 


In a plate problem inplane strains (e^, Y^y) and thus 
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stresses (a^^, cTy, are the functions of the height, while out of 
plane strains Yy^) and stresses (t^^, Ty^) are constant across the 
thickness. And the stress strain relationship, 

{a} = [DD]{e} (2.10a) 


[DD] 


1 V 0 0 0 

V 1 0 0 0 



1-v 

2 

0 


0 

1-v 

2k 


0 

0 


0 0 0 


1 - 

2k 


{2.10b) 


Where 

i<y} = ^xy > {e} 


In eq(2.10b) k=1.2 is a factor for improving the 

performance of the analysis, and this accounts for non uniform 
shear stress distribution. The stress and strain vectors have 
five above mentioned components. 

b. Yield Criteria 

In the present analysis Prandtl-Reuss relation is made use 
to determine the onset of yield. 
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( 2 . 11 ) 


^lol + ol - a^Oy + 3x^ +3 tL + 

c. Fraction Of Plasticity 

As the load is increased, it is known, that the outer 
most fibres become plastic first and the plasticity spreads in 
to the core and across the plate. At the same time whole 
thickness is not plastic. So to find above what height the plate 
goes plastic, the height of elasticity is calculated as (fig 2.1) , 

2 2 

^ ^ yield ~ 3 (Xxz Jl (2.12) 

ol + ol - o^Oy + 3-cly^ 2 


If the numerator of the term inside the root sign goes 
negative, the whole thickness goes plastic. 

d. Flow Rule 

Plastic materials are governed by flow rule between stress 
and plastic strain increment. The flow rule is obtained by the 
normality principle. The relation between effective stress and 
strain in simple terms can be understood by the power law in one 
direction. The power law is an experimental agreement of the 
material properties in the plastic range. 

a = (2.13) 


Where K and n are material constant and calculated 
experimentally. The eg (2. 13) governs the relation only after 
yield point. The power law is shown in fig (2. 3) . 
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e. Plastic Stress -Strain Matrix ( [DD]*** ) 

If a point in the plate goes plastic then the plastic 
stress strain matrix is used at that point to calculate the 
incremental stresses from incremental strains in the plastic 
region. 


{6a} = [DDl^p{be.} 


(2.14) 


Where 


[DD] = [DD] - { 2 } 

AA + {a}'^[DD] {a} 


Here [DD] is the property matrix, AA. is the slope of stress- 
strain curve governed by elastoplstic strain stress relation and 
{a} is the flow vector which depends on the plasticity criteria. 

f. Iterative scheme 

The unbalance forces to be calculated by commonly used 
equation, 


{f} = ’[BS]^{a}ciF 


(2.15) 


does not work in this mixed formulation because of the 
following reasons . 

1. The eq(2.i5) can be used if the following expression 
satisfies the stiffness matrix of an element. 
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[ic]® = 


[BB]'^[DD] [BB]dV 


(2.16) 


2. Also if the above equation is assembled to form global 
stiffness matrix for this particular robust element, the global 
stiffness matrix goes singular. 

Answer to these difficulties is found in the mixed 
formulation which incorporates shear forces also as degree of 
freedom. Inclusion of such degree of freedom avoids the defects 
of singularity and locking. Presence of changes the character 
of stiffness matrix, while the same is not present in 
eq(2.16).The eq(2.7) is derived from "mixed formulation" while 
eq(2.15) and eq(2.16) represent the "irreducible formulation". 

So leaving eq(2.15) and eq(2.l6) aside, another iterative 
scheme is tried which make use of original elemental stiffness 
matrix. 

For an element, if it goes plastic at a given load then the 
incremental deformations should satisfy the following equation. 

[iC]®^{5u} = {6f} (2.17) 

Here [K]®p is the elastoplastic stiffness matrix. Since this 
matrix is a function of [DD] matrix and thus a function of 
stresses, the global [K] matrix needs updating. So to avoid such 
recalculations use of the iterative scheme of eq(2.l7) is 
implemented in the following form. 
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= {5/} +[[i?]®- [J?]®^]{6u} 


Now displacement vector on left hand side is treated as 
displacement vector obtained after n*^ iteration while the 
displacement vector on right hand side is the vector before the 
iteration and is used to calculate unbalance forces. Now for 
showing the iterative scheme the eq(2.18) can be rewritten as, 

[iC]^{6u}„= {6^} (2.18) 


The geometrical representation of eq(2.18b) is shown in the 
fig(2 .4) . 

g. Convergence Criterion 

The "error" in each iteration is represented as follows, 

e = {uj - 

The norm of this error vector will be, 

iej = [e^’e] ^ 


And convergence is said to be achieved when. 




(2.19) 


Where a is the convergence limit. 
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h. Effective Plastic Stiffness [K]®p Matrix 

[K] matrix is the plastic stiffness matrix for an 
element. Since the whole plate may not go plastic and some 
portion may still remain elastic, the calculation of this matrix 
have to made carefully. 

Eq(2.14) shows incremental stress strain [DD]®p matrix for 
full plasticity and eq(2.10b) shows stress strain matrix [DD] for 
full elasticity. In [K] ® matrix eqn(2.8) the material properties 
are used at two places in form of [D] and [C^] matrices. Now 
matrices [DD] and [DD]®p are defined in the following way. 


[DD] 


[DD^] 0 

0 [DD^] 


■[DDJ 

[DD] = 

0 


0 ' 
[DDJ 


Where [DD,] , [DD,]®p are 3x3 and [DDj] , [DD 2 ]®p are 2x2 matrices. 

Now effective plastic stress strain matrix can be written in the 
following way. 


* [jj - ( 2 . 20 ) 

Which represents the gauss numerical integration with 7 
gauss points of the element. [DD,]®*^ and hp are also different for 
different gauss points. 

And 

- (h-2h^)[DD,]^P]w, (2 2 ^) 
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Fitting [D]®p and [C,]®p in place of [D] and [C,] in the 
eq(2.6), The [K]®^ matrix in place of [K]® matrix is obtained. 

The loading and unloading parts change the characteristic 
of stress strain curve and the yield stress. This aspect is also 
included and discussed in next Chapter with other aspects in 
implementation of the computer code. 



CHAi'riSK 3 


IMPLEMENTATION OF FINITE ELEMENT CODE 

This chapter explains the salient aspects of numerical 
method in developing the FEM code, including element selection. 


3.1 Element Selection 

Triangular elements are easier for discretizing complex 
geometries. In literature for different types of problems 
different type of elements are suggested as they should be free 
from defects of locking and singularity. Also the performance 
of the element under different type of loadings and boundary 
conditions should be stable. "Robust triangular element" with 
suitable nodal variables satisfy the stability conditions. 

The choice of the shape functions and variables should be 
such that the Babuska-Brezzi [Zienkiewicz and Lefebvre, 1988] 
conditions are satisfied by the system and thus, stability is 
ensured. The most vital and strictly necessary condition for the 
non- singularity of the system given by eqn(2.8) is that, 

(3.1a) 


Ct = 


24 


(3.1b) 


P 



n 


W 


Where ng,ng and stand for the number of variables for 
rotation, shear force and lateral displacements respectively in 
an element. If either of the two conditions are violated, the 
locking behaviour will occur in the numerical solution. 

For the dynamic analysis if an implicit time integration 
scheme is used, it is necessary to use higher order finite 
elements and a consistent mass matrix. The higher order elements 
are effective in the representation of bending behaviour. These 
elements should be employed with consistent load vector, in that 
the midside and corner nodes are subjected to their appropriate 
load contribution in the analysis. 

In the mixed finite element formulation, triangular element 
(T6/3) has six boundary nodes which are used to determine a 
quadratic variation of 0 and w [Zienkiewicz and Lefebvre, 1988] . 
Shear forces are determined by a linear field defined by three 
internal nodes . 

When standard patch tests [Zienkiewicz and 
Lefebvre, 1988 /Zienkiewicz et al.,19863 are used with this 
element, the results indicate that, the element has incipient 
locking possibility and thus it is not a satisfactory element. 
But the examination of results indicate that the satisfaction of 
the requirements can be achieved by the addition of the 6 
variables in the interior of the element. The three internal 


0,vj 2,1 d, o. f 


s 2 d. o, f 
s,0 2,2 d c 

O C; - 'O 

(O) T6/3 '£4 3 f , 



O ° 'C ,. 

□ □ 

O o ^-) 

(b) T6/3 B3 (30 d, o, f.) 
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nodes created have the shape functions used of the form L^^L 2 L 3 , 
provide a new element (T6/3 B3) which passed the same patch test 
successfully in all cases [Zienkiewicz and Lefebvre, 1988] . 

This particular element is chosen for the present analysis. 
It has 30 degree of freedom, and after the condensation of shear 
force terms the number of variables reduces to 24 . 


3.2 Time Integration scheme 

The eqn(2.8) is integrated by a numerical procedure. The 
integration is based on two ideas. 

1. Static equilibrium including inertia force is achieved at 
discrete time intervals At. 

(3.2) 


2. Variation of displacement, velocity and acceleration within 
each time step is assumed. This assumption decides the accuracy 
and stability of the solution procedure. 

The idea no. 2 decides the time integration scheme. The 
choice of a scheme depends on the finite element idealization. 
However, the finite element idealization to be chosen depends in 
turn on the actual wave velocity of the medium to be analyzed. 
It follows, therefore the selection of an appropriate finite 
element idealization of a problem and the choice of effective 
integration scheme for the response solution are closely related 
and must be considered together . 

Out of many schemes "Newmark's implicit scheme, is most 
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accurate and unconditionally stable. 


3.2a Newmark Implicit Scheme 

The Newmark 's constant acceleration scheme is written as 
follows [Bathe, 1990] , 

= <^u + [(1-5) At (3.3a) 

t*^tu = + [ (A-ct) ^f.2 


Where a and 6 are parameters that control integration 
accuracy and stability. Newmark originally proposed an 
unconditional stable scheme, in which case 6 = 0.5 and 
a = 0.25. 

Using eqn(3.2) and eqn(3.3) final form can be obtained as, 

+ [M] {ao^U + + a^^'u) (3.4a) 


-^U) - - aj^U 


(3.4b) 


t+Aty 


(3.4c) 


Where , 


I 
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The implicit scheme, as 
formulation requires inversion 
purpose Cholesky decomposition 
efficient . 



can be seen from the final 
of [K+apM] matrix. For this 
technique is computationally 


3.2b Time step 

After deciding on time integration scheme, it is necessary 
to use a time step suitable with finite element mesh. If the 
critical wavelength in the medium is denoted by L^, then the time 
step can be decided as, 

At = 

c n 

And the effective length of a finite element should be 

Lg = c At 


Where c is wave speed and 'n' is the number of nodes per 
critical wavelength. 

The effective wavelength and corresponding time step must 
be able to represent the complete bending wave accurately and is 
chosen appropriately depending on the kind of element 
idealization and time integration scheme used. It is reported 
that at least 8 nodes per shortest wavelength are required in 
order to produce all artifacts of wave propagation [Seron,1990] . 
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However, for higher order (parabolic and cubic) continuum 
elements the time step has to be further reduced, because the 
interior nodes are stiffen than the corner nodes [Bathe ,1990]. 
In addition for T6/3 B3 , the internal 0 variables are quartic 
functions. Also in the plate elements the flexural modes also 
affect the time step. In the actual computation appropriate time 
step is to be selected. 


3.3 Implementation Of Plastic And Dynamic Code 

At any time step the governing equation of the dynamic 
problem can be defined as eqn(2.8) . Mass matrix [M] does not get 
affected by the plasticity. So for art incremental load the 
equation can be written in the elastoplastic region. 

[MlHbU} + = {6f} (3.5a) 


Again rewriting the above equation as we did for eqn{2.17) , 
the equation acquires the following form, 


[Af3®{SU} + 


(3.5b) 


Right hand side shows the sum of the original load and the 
unbalance. The unbalance arises because of the difference in 
elastic and plastic stress strain matrices. For every time step 
this equation is- balanced to get an equilibrium displacement 
satisfying the plasticity in the element. Initial unbalance is 
taken to be zero and then a trial displacement vector is 
calculated. This process is continued until both side equalize. 


30 



The advantage of the eqn(3.5) is that the global stiffness 
matrix is not calculated repeatedly, and thus saving computer 
time . 


3.4 Loading And Unloading 

If a material becomes plastic while loading and then 
unloaded, it behaves purely elastically. To implement this aspect 
the history of the load is to be traced to keep track of the 
strains in the material. Once the material is unloaded the yield 
point also changes. So in case of plastic loading and unloading 
the yield point of every gauss point is different. In case of 
dynamic loading the arrival of plasticity may have two effects 
on the stability. 

1. The displacements are higher than those in elastic case and 
thus functional aspect of the design gets affected. 

2. The plasticity absorbs much energy, thus the impact of 
dynamic loading may reduce. Also in case of loading and unloading 
the energy which comes under hysteresis loop gets released, so 
the damping and stability in the system may increase. 


3.5 Con^uter Implementation 

The computer code developed by Dindore [1994] has been 
modified to take the elastoplastic nature of the material in the 
present work. 

The present work deals with two type of loadings, uniform 
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distributed loading and moment loading on the boundary. Plastic 
code first determines the stresses and if these are still in the 
elastic region, the plastic part is bypassed. Also if ultimate 
stress is reached, the programme stops. If the load is in plastic 
region the load above the yield point is divided into 20 equal 
part and then applied in small steps. The core size of elasticity 
(hp) and stresses are updated after each force iteration to save 
much computer time. The maximum error due to such scheme can be 
justified in following way. 

Maximum stress generated < a ,, < 1.5 o ■ 

ult yield 

so, stress increment < . °y. = o. 025 o_j. 


So the percentage of elastic height that will be missed due 
to above mentioned assumption in each force iteration is, 


bh^ < 


0.0 2 5 03, 


xlOO =2.5% 


This error gets compensated to the extent because of 
following reasons, 

1 . Out of plane stresses are not varying with height . 

2 . Weightages of various gauss points average out the 
error. 

3. Averaging out if the number of elements are more. 

Out of seven gauss points any one or more may be elastic 
or plastic upto any height in the thickness. The results are 
shown as contour plots to show extent of plasticity on the plate 
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CHAPli^ 4 


RESULTS AND DISCUSSIONS 


This Chapter is divided into two parts. Former consists of 
validation of computer code while later analyzes the 
interdependence of the plastic and dynamic behaviour. 


4.1 Validation of code 

Validation of computer code is done on following cases, 

case 1. Stresses and deflections for the purely elastic problem 
are calculated for standard cases (Timeshenko, 1959) . A 
satisfactorily agreement (within 1.5 %) is observed with 
theoretical results. Mesh is used as in fig (4. 2) . 

case 2. For qualitative validation, the plate is subjected to 
different loads and boundary conditions. The spread of 
plastic zone was observed and the same was more in the 
cases where either geometric curvature was more or there 
was a discontinuity in boundary conditions . 

case 3. For quantitative validation, a standard problem is taken 
up (Hinton, 1984) . In this a plate fixed in all edges, 
subjected to a uniform load is considered (fig 4.1) . The 
material and geometric properties are as follows. 
Modulus of elasticity, E = 30000 MPa 

Poisson's ratio, v = 0.30 

Density P — 0.003 x 10^ Kg/m^ 
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Thickness of plate, h = 0.20 m 

Side of the square plate, L = 6 . 0 m 

slope in plastic region, K = 0.01 x E 

Yield stress, 30.0 MPa 

Uniformly distributed load for q = 0.30, 0.35 ,0.40 MN/m^ 

The present finite element analysis consists of dividing 
the square plate by a uniform mesh as shown in fig (4. 2) . As can 
be seen quarter plate is also used by applying symmetric 
conditions at the boundaries. 

Fig (4. 3) presents the deflection at the centre of the plate 
as a function of load. The figure shows the results of 
considering the mesh as a full plate and as quarter plate (which 
is equivalent to a finer mesh with elements of half size) 
together with those semi-loof analysis [Hinton and Owen, 1980] 
and Hinton and Owen, 1984. It can be seen, that the present finite 
element results are in close agreement with semi-loof analysis. 
Fig(4.4-4.5) show the spread of plastic zone as contours plots 
at various load levels. At each load level the different contours 
correspond to percentage of plastic height. These plots also 
agree well with the results of Hinton (fig 4.4a) except toward 
peak load (q=0.40) Hinton's results behave more plastically. This 
happens because of layered approach used by Hinton. An inherent 
error of 12.5 percent is induced in dividing the thickness into 
8 layers. And at peak loads the present results agree with those 
of Hinton's in this error range. 

Thus, it can be seen that the results by the present code 
agrees well with the available quasi-static elastoplastic 

results . 



4,2 Pldstic d.xid dynamic bcliavioim 


This section presents elastoplastic dynamic results 
obtained by the present approach. The variation of dynamic load 
is assumed as follows, 


Q = 



cos 


n t" 
RT 


t ^ RT 


(4.1) 


For load # 1, 


Q = Qo t > RT 


And for load # 2, 


Q = 


Qo cos 


' Tt ( t - RT) 
2 RT 


t > RT 


Where RT is the rising time for the load to reach the peak. 
The forcing functions are shown in fig(4.6). The element size 
and time steps are choosen in such a way so that the solution 
remains stable. 

To study the interdependence of plastic and dynamic 
behaviour, the analysis is divided into three parts. 

I. The amount of plasticity which is the ratio of effective 
stress to yield stress in the most plastic gauss point, is shown 
against time for various rise time (fig 4.7) . The load no. 1 is 
applied for this case. This is compared with the amount of 
plasticity in static case and also with onset of plasticity in 
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that elainsnt . For this case the amount of plasticity increases 
as the rise time is decreased while keeping other parameters 
same. By reducing the rise time the inertial forces play a 
significant role in inducing more plasticity. 

II . For a given rise time the dynamic response is analyzed for 
purely elastic and elastoplastic cases. The load no. 2 is applied 
in this case. The deflection is also seen for different loads 
with same rising time. The deflection at the centre is plotted 
against time steps. Initial maximum deflation is more than the 
purely elastic case. But at latter time the delection is less 
(fig 4.8). For this case the dynamic response in elastoplastic 
case is more damped than the purely elastic case. The response 
becomes more dampened as the plasticity increases in the plate. 
The dissipation of energy in form of hysteresis losses accounts 
for this damping. 

III. The number of gauss points which have gone plastic, are 
plotted against the time for different loads and given rise time. 
The load no. 2 is applied in this case (fig 4.9). For this case 
because of strain hardening the number of elastic gauss points 
is increasing during unloading. And at the next peak load the 
number of plastic gauss points is reducing than the previous one. 



-q 



rig (4,1) 


Where q is the uniformly distributed load. 


Fig( 4.1 ) 
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Fig (4. 7a): Effect of dynamic loading on the plasticity 
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Fig(4.7t)): Effect of dynamic loading on plasticity 
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CHAPTj^ 5 


CONCLUSIONS AND SCOPE FOR FUTURE WORK 

5 . 1 Conclusions 


Based on the present finite elastoplastic and dynamic 
analysis of a plate under dynamic and static loading, following 
conclusions can be verified, 

1 . In case of dynamic loading on the plate the plasticity is 
more than it was in purely static load with the same peak value. 

2. The plasticity keeps on increasing as the rise time is 
decreased. And for a sufficient longer rise time there is no 
additional increase in the amount of plasticity over the static 
case . 


3 . As the plasticity is induced in the system, the energy is 
dissipated which acts as a dampener on the dynamic response and 
thus the oscillations reduces or vanishes. 


4. The deflection during first peak load is more in 
elastoplastic case than purely elastic case. And in latter load 
peaks because of damping the deflection comes down for 
elstoplastic case. 


5 . 2 Future Scope Of Work 


The mesh can be made finer and the time step can be made 
smaller to observe the sharper impact loading and plasticity more 

accurately. CENTRSL LiBRAOT 

I » ' AAMPyR 


A 1 Z I tl* 5; 
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Several other loadings and different geometrical 
configurations can be tried. Also the cases of cyclic loading can 
be studied to see the permanent deflection after a given number 
of cycles. 

Normally the small cracks are always present in the 
material, but that can not be modelled in the finite element. So 
to model the cracks we can make use of statistical method with 
varying yield stresses at different gauss points. 

The effect of plasticity acts as an dampener in the dynamic 
response. The study can be useful where the system vibrations are 
to be avoided at the cost of permanent damage of some parts of 
the system. 

For dynamic loading the time step can be selected in a 
optimum manner so that the computational time is less with 
considerable accuracy and stability. 
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AJESTTLACT 


In practical cases when a pLate is subjected to a dynamic 
loading , the geometric artd nrateirial linearity may not be 
preserved. This study is aiied to investigate the effect of 
plasticity in the dynamic response . _As there are no closed form 
solutions available owing to the ccomplexity of the problem 
, numerical methods are used to s ol ve the problem. 

For this computational work an efficient FEM code is 
developed using the "robust trdang-ul ar element" developed by 
Zienkiewiz [1991] . "Mindlin pLate ttieory" in the form of mixed 
formulation is used for tnodeblcng pLate bending. The inertia 
effect of the element is repie s&ntedl by consistent mass matrix. 

As the stresses are var^^img along the thickness, the 
plasticity should be applied togiether with elasticity for every 
iteration. This introduces tbe concept of "elastic height" , which 
decides the stresses In the plates. And to find out the 
elastoplastic stiffness matrices., th_e integration should be done 
in both elastic core and. plastic outer region. 

The computer code is firs t ms edL to solve a standard problem 
and after validation it was ms ed boa: dynamic cases of impact 
loading. The dynamic loading increas^es the plasticity and in turn 
the plasticity in the material acts as an dampener to the dynamic 
oscillatory response . The interdleperndence of plastic and dynamic 
behaviours is analyzed in tine present work. 



1.1 GENERAL INTRODUCTION 


Plasticity in the plate is different from a simple two 
dimension analysis of the plasticity. The main difference is the 
variation of stress along the thickness. As the first 
approximation of the plate stress and strain in the thickness 
direction are neglected [Reissner Mindlin plate f oirmulation] . And 
if the plate is thin then the out of plane stresses are 

also neglected. For a thick plate in-plane-stresses cTy, 
varies with the height while out of plane stresses remain almost 
constant along the height. This variation of the stresses affects 
the size of the elastic region in the core of the plate. 

In general plate formulation is made in terms of lateral 
forces and moments. In plastic case since the in-plane stresses 
are different in the elastic core and outer plastic zone, the 
calculation of the stress resultant moments becomes a function 
of the size of elastic core. The whole thickness is not plastic, 
proper weightages are to be given to the elastic and plastic 
stiffness matrices to calculate effective stiffness matrix. 
Calculations of unbalances after the integration of dynamic 
analysis and the plastic analysis is to be done to incorporate 
both the aspects properly. These taken together increases the 
complexity of the problem. 


1.2 LITERATURE SURVEY 


The thin plate theory developed by Kirchhoff is modified for 
thick plates by Reissner and Mindlin by relaxing the assumption 
that the normals to the midplane remain normal. This was done 
since the shear force terms can not be neglected in many loading 
and support conditions. Reissner and Mindlin theories can be 
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extended to thick plates and thus known as Reissner-Mindlin 
formulation. Regarding Finite Element Method a difference 
between thin plate and thick plate analysis is that in former it 
is possible to represent the state of deformation by one quantity 
deflection w and thus C’ continuity of w is needed, while in 
latter w can have C° continuity along with rotation and shear 
degree of freedom. 

While analytical solutions remain difficult for thick 
plates, the advent of Finite Element Method made thick plate 
theory simpler to implement. Thick plate theory on the imposition 
of additional assumptions behaves as thin plate theory. It is 
more natural to start with thick plate fomtiulation and make 
transition to thin plate theory by imposition of certain relation 
between degrees of freedom [Zienkiewicz and Taylor, 1991] . In 
thick plate formulation instead of one approximate function for 
w ( for calculation of rotational angle ) , independent 
interpolation functions for deflection w and rotation 6 are 
taken. During formulation, if the shear forces are eliminated, 
the formulation is known as "irreducible formulation" and when 
functions for shear forces are also retained the formulation is 
known as "hybrid formulation" . The latter formulation gives more 
degree of freedom and so it is preferred to work with in most of 
the cases [Zienkiewicz and Taylor, 1991] 

Another aspect is the selection of the element to perform 
satisfactorily for thick plates as well as thin plates. Many 
elements are attempted to avoid the problems of locking and 
singularity but the results of "triangular robust element" are 
found favourable in this regard [Zienkiewicz and Lef ebvre, 1988] . 

For determination of plasticity three things are required: 
an yield criteria, flow rule and plastic stress strain law. A 
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fundamental difference between plastic and elastic analysis is 
that in elastic solution total stress can be evaluated from the 
total strain alone, while in the plastic response evaluation of 
stress strain at any time depends on the stress strain history. 
This depends on the experimental results which are termed as 
yield criteria. Commonly used yield criterion are continuous one 
but piecewise linear yield criterion can also be used [Giulio 
Maier et al.,1987]. 

Yield criterion defines the onset of plasticity and flow 
rule governs the relationship between stress and plastic strain 
increment after the yield point. Flow rule requires the 
definition of a potential which is taken to be the same as yield 
surface for the case of associated plasticity. The plastic stress 
strain law has been defined for nonlinear, linear and perfectly 
plastic cases [Yamada et al., 1968]. In latter two cases the 
total stress strain curve is bilinear and in former case usually 
some power law satisfying the experimental results is defined. 
Flow rule, thus can be used in plastic cases by forming a plastic 
potential and then defining the plastic strains by the same. When 
the plastic potential surface is the same as the yield surface, 
the plasticity is termed as "associated" otherwise the plasticity 
remains "unassociated" . The flow rules for loading and unloading 
are not same, and that makes the problem difficult. 

Out of many iterative schemes to solve the above non-linear 
problem, Newton-Raphson method is found favourable in terms of 
ease in mathematical formulation and in computational effort 
[Simo and Taylor, 1986] . In Newton-Raphson method, incremental 
strain is multiplied to plastic stress strain matrix to get the 
stresses, while in Runge-Kutta method the incremental strain is 
subdivided in many subincrements then individual subincremented 
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strains are multiplied to their corresponding stress strain 
matrices to get subincremented stress. These stresses are added 
up to get incremental stress . This is known as 
"subincrementation" [Nayak,1971; Bathe, 1990] . 

For plasticity in the Mindlin plates, it is suggested to 
use full three dimensional analysis of stress field. Inclusion 
of shear terms in the yield function allows the spread of 
plasticity from the extreme fibre to entire plate thickness 
[Hinton and Owen, 1980] . A very thick plate can be divided in many 
layers of thin plates and with the Mindlin plate theory, the 
problem is solved [Hinton and Owen, 1980]. But the triangular 
robust element is proved favourable for transition from thick to 
thin plates, and the non layered approach can be used for 
moderately thick plates [Zienkiewicz and Lefebvre, 1988] . 
Elastoplastic numerical experiments has been carried out adopting 
a full three dimensional version and different criteria of yield 
in the past [Hinton and Owen, 1984] .These results show excellent 
agreement in static elastic analysis and the prediction of 
vibrational frequencies and buckling loads. These results forms 
the basis of comparison for the present work. 

In elastodynamic, a detailed analyses of various time 
integration schemes from point of view of accuracy and speed, 
were reported by [Seron et al.,1990; Zienkiewicz and Taylor, 1991; 
Bathe, 1990] . It is concluded that, even though central difference 
scheme is the best in terms of computational speed and cost, 
Newmark's constant average scheme gives better accuracy and 
offers unconditional stability. Whereas, the time step size has 
negligible effects when central difference scheme is employed 
[Bazant and Celep,1983]. 
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In elastoplastic dynamic analysis based on Mindlin plate, 
the work started quite late. Most of the work is in the area of 
concrete structure under dynamic loading [Kujawski et al.,1984] . 

Experimental information on material behaviour is scarce 
for most structural materials [Hinton and Owen, 1980] . 
Instantaneous yield stress is significantly influenced by the 
rate of straining, also the value of elasticity modulus is found 
to be dependent on the strain rate for structural material those 
varies from the materials of limited ductility to brittle elastic 
behaviours. For this purpose a better understanding of observed 
phenomenon and underlying microscopic behaviour is required 
[Hinton and Owen, 1980] . 

In loading and unloading the yield point of the material 
changes because of strain hardening. A model for yield limit 
degradation has been developed by Bazant [1987] . Other study 
[Panzeca et al.,1987] deals with a transformation matrix that 
takes care of yield point variation according to the loading and 
history. 

Application of numerical model to dynamic flow of steel 
specimen has been done by Casadei et al , [1987] . The results 
compared favourably with the experimental result. For beams and 
rods a code validation case study is done by Marques [1987] . An 
elasto-plastic one degree of freedom numerical model was 
presented for steel slab under blast loading [Mazza et al.,1987] . 
The results were compared with static elastoplastic loading. Also 
both of the results were compared with experimental data. It was 
found out that in the dynamic loading conditions the 
displacements were intensified at the plastic zone compared with 
those of static loading. In the area of explosive welding also 
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the study was carried out [Ante, 1987] . In this case welding of 
two concentric tubes were considered while explosive was put 
inside the inner tube. The specimen experienced both mechanical 
as well as thermal shock. In this problem plasticity was induced 
by thermal load also. 


PRESENT WORK 

The present work deals with the study of the effect of 
plasticity in plates under various type of loadings both static 
and dynamic . 

A mixed finite element analysis using Reissner Mindlin 
plate theory is used to model the plate bending. To avoid locking 
and singularity a triangular robust element (T6/3 B3) [Zienkiewicz 
and Taylor, 1991] is used. 

Chapter 2 deals with the theoretical part of the 
formulation of plate bending and plasticity. Chapter 3 deals with 
the finite element implementation of the plastic analysis and its 
integration with dynamic analysis. Chapter 4 describes the 
validation of the code and the behaviour of plasticity in dynamic 
loading and effect of plasticity on the dynamic response. Chapter 
5 presents the conclusions and the suggestions for future work. 
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CHAPTER 2 


BASIC EQUATIONS OF ELASTOPLASTIC DYNAMIC PLATE ANALYSIS 

2.1 INTRODUCTION 


This chapter presents the basic equations governing the 
plate under elastoplastic and dynamic conditions. As explained 
in earlier chapter, by appropriate element selection and 
formulation Mindlin plate theory can be used for thick and thin 
plates . 

The Mindlin plate theory makes the following assumptions, 

1. The strain and the stress components in Z-direction are 
negligible ( cr^ and «= 0 ) (fig 2.1), where Z-direction 
represents the thickness direction. 

2. The normals to the middle plane of the plate remains straight 
during deformation. 

Applying above assumptions, static equilibrium equations 
for plate bending problem can be obtained [Timoshenko, 1959] in 
terms of midplane rotations {6j(,6y) , lateral displacements (w) of 
the midplane and the shear forces (S^^, S^) at the cross-section. 
Sign convention and the variables are shown in fig{2.1). 

Stress resultant moments and shear forces are defined in 
terms of stress variables as follows: 
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(2.1a) 


M = 


/. 


h/2 

h/2 


z dz 


xy 


(2.1b) 




= dz 


rh/2 


(2.1c) 


For the dynamic analysis inertia forces due to mass and 
moment of inertia are considered in equilibrium equations. 
Formulating equilibrium equation using d' Alembert's principle, 
the following equations are obtained. 


dx 




(2.2a) 


dM 

dx 




(2.2b) 


+ Q - phw = 0 (2.2c) 

dx oy 


Strains at a height of z from the centre plane are given as, 
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(2.3a) 


36 . 

^ ^ dx 




dy dx _ 


a , dw 

" to 


Yyz = % + 


8w 

dy 


(2.3b) 

(2.3c) 


2.2 THE MIXED FINITE ELEMENT FORMULATION 

The finite element formulation involves three sets of 
variables 0,S,w. These variables at any point within the element 
are expressed in terms of nodal parameters (6^^, 6y, Sy,w) of the 
element using shape functions. 



0 = 




NqQ 


(2.4a) 



(2.4b) 


w = (2 .4c) 

Where 0,S,w are the nodal values of respective variables. 
Now using eg(2.3) and eq(2.l) the eg(2.2) can be rewritten 


1 , 
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as follows. 



(2.5a) 


L^DLQ - S -p— 6 = 0 

r - ^ 


e - q5 + vw = 0 (2.5b) 


+ g - phw =0 (2 . 5c) 


Where L,V operators and matrices are given by, 


dx 


L = 


0 


d 

dy 


0 


ai 

a 


dx 

dy 

/ V = 

a 

a 


. 3y. 

dx 




^1 


j\i 0 

Gh 0 1 


D 


fl V 


Eh^ 

12 (l-v2) 


V 

0 


1 

0 


0 

0 

1-v 

2 


Now by Galerkin method, using eq(2.5) and eq(2.4), and 
then integrating over the area of the element, the following 
equations are obtained. 


A0 + B5 + £8 = fg 


(2.6a) 


B^d - PS + Cw = 0 


(2.6b) 
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(2.6c) 


C^S + Fw - 


Where , 


A = f^(LNQ)^D{LNQ)dQ , B = f^Ne^N^dd 


C = fj^JivNjdCl , P = j^Njc^N^dQ 


fQ = f Ne" T dr , f„= ( nJ g dOl 

jTj- jq 


E = Nq dCl , F = p hJ^Nj dQ 


In the above equations the vector T represents two 
moment components imposed on the boundary by the traction and q 
is the intensity of the imposed lateral force. The shape 
functions Ng,Ny are chosen to have Cp continuity, but because of 
shear forces are defined only for inside nodes of the element, 
R can be discontinuous between the elements. Such a choice 

S 

allows S to be defined locally for each element and thus be 
eliminated when P is a non- singular matrix. Thus the problem 
includes only 6 and w as variables . 

The final formulation in terms of 6 and w, after 
condensing S, becomes : 
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{ 2 . 1 ) 


A 4 BP-^C 

e 

+ 

E 

o’ 

e 




w 


0 

F 

w 




In short eg (2. 7) can be written as, 


[ir]®{u} 4 M®{u} = {F}® 


( 2 . 8 ) 


2.3 PLASTIC ANALYSIS 


When the stresses go beyond the yield value, then stiffness 
matrix changes depending on the individual stress component, 
effective stress and flow rule. In such cases no direct solution 
can be used, so iterative schemes are used to solve the problem. 

To do the non-linear problems, the following steps are 

used: 

a. Calculation of stress and strain 

After obtaining the maximum possible elastic solution, the 
correspondig strains are calculated by the following formula 
which is obtained by substituting eg (2. 3) using the variables 
from eg (2 .4) . 


{e} 


z{LNq) 


0 

viV, 


{u} = [BB] {u} 


(2.9) 


In a plate problem inplane strains y^y) and thus 
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stresses (cT^, CTy, are the functions of the height, while out of 
plane strains Yy^) and stresses (t^^, are constant across the 
thickness. And the stress strain relationship, 

{a} = [DD] ie} {2.10a) 

1 V 0 0 

V 1 0 0 

0 0 0 
2 

0 0 0 — 

2 . 

0 0 0 0 




Where 



In eg (2.10b) k=1.2 is a factor for improving the 

performance of the analysis, and this accounts for non uniform 
shear stress distribution. The stress and strain vectors have 
five above mentioned components. 

b. Yield Criteria 

In the present analysis Prandtl-Reuss relation is made use 
to determine the onset of yield. 
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^ol + al - o^Oy + 3xly +3xL + 3^^ = 


(2.11) 


c. Fraction Of Plasticity 

As the load is increased, it is known, that the outer 
most fibres become plastic first and the plasticity spreads in 
to the core and across the plate. At the same time whole 
thickness is not plastic. So to find above what height the plate 
goes plastic, the height of elasticity is calculated as (fig 2.1) , 


= 


^ yield ^ '^yz^ 




O^Oy + 3'lxy 


1/2 


h 

2 


( 2 . 12 : 


If the numerator of the term inside the root sign goes 
negative, the whole thickness goes plastic. 

d. Flow Rule 

Plastic materials are governed by flow rule between stress 
and plastic strain increment. The flow rule is obtained by the 
normality principle. The relation between effective stress and 
strain in simple terms can be understood by the power law in one 
direction. The power law is an experimental agreement of the 
material properties in the plastic range. 

a = Ke^ (2.13) 


Where K and n are material constant and calculated 
experimentally. The eq(2.13) governs the relation only after 
yield point. The power law is shown in fig (2. 3) . 
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e. Plastic Stress-Strain Matrix ( [DD]*p ) 

If a point in the plate goes plastic then the plastic 
stress strain matrix is used at that point to calculate the 
incremental stresses from incremental strains in the plastic 
region. 


{5a} = [DD\^P{be} 


(2.14) 


Where 


[DD]^^ = [DD] - t-P-P] ^3} {a}5_[g.D]_ 
A4 + {aVlDD] ia] 


Here [DD] is the property matrix, AA is the slope of stress- 
strain curve governed by elastoplstic strain stress relation and 
{a} is the flow vector which depends on the plasticity criteria. 

f. Iterative scheme 

The unbalance forces to be calculated by commonly used 
equation, ' 

{f} = (2.15) 

does not work in this mixed formulation because of the 
following reasons. 

1. The eq{2.15) can be used if the following expression 
satisfies the stiffness matrix of an element. 
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[iC] ® = J [BB] "^[DD] [BBl dV 


(2.16) 


2 . Also if the above equation is assembled to form global 
stiffness matrix for this particular robust element, the global 
stiffness matrix goes singular. 

Answer to these difficulties is found in the mixed 
formulation which incorporates shear forces also as degree of 
freedom. Inclusion of such degree of freedom avoids the defects 
of singularity and locking. Presence of Ng changes the character 
of stiffness matrix, while the same is not present in 
eq(2.16).The eq(2.7) is derived from "mixed formulation" while 
eq(2.15) and eq{2.l6) represent the "irreducible f oinmulation" . 

So leaving eq(2.15) and eq(2.l6) aside, another iterative 
scheme is tried which make use of original elemental stiffness 
matrix. 

For an element, if it goes plastic at a given load then the 
incremental deformations should satisfy the following equation. 

= {6f} (2.17) 

Here [K]®p is the elastoplastic stiffness matrix. Since this 
matrix is a function of [DD] matrix and thus a function of 
stresses, the global [K] matrix needs updating. So to avoid such 
recalculations use of the iterative scheme of eq(2.17) is 
implemented in the following form. 
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[m^ibu] = { 5 /} +[m^- 


Now displacement vector on left hand side is treated as 
displacement vector obtained after n^^ iteration while the 
displacement vector on right hand side is the vector before the 
iteration and is used to calculate unbalance forces. Now for 
showing the iterative scheme the eq(2.18) can be rewritten as, 

[iC]^{6u}„= {6/} (2.18) 


The geometrical representation of eq(2.18b) is shown in the 
fig(2 .4) . 

g. Convergence Criterion 

The "error" in each iteration is represented as follows, 

e = {uj - 

The norm of this error vector will be, 

J. 

|e| = [e ^e] ^ 


And convergence is said to be achieved when. 



(2.19) 


Where a is the convergence limit. 
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h. Effective Plastic Stiffness [K]*p Matrix 

[K]®P matrix is the plastic stiffness matrix for an 
element. Since the whole plate may not go plastic and some 
portion may still remain elastic, the calculation of this matrix 
have to made carefully. 

Eq(2.14) shows incremental stress strain [DD]®^ matrix for 
full plasticity and eg(2.10b) shows stress strain matrix [DD] for 
full elasticity. In [K]® matrix eqn(2.8) the material properties 
are used at two places in form of [D] and [C^] matrices. Now 
matrices [DD] and [DD]®p are defined in the following way. 


[DD] 


[DDJ 0 
0 [DD^] 


[DD] 


[DDJ 0 

0 [DD^] 


Where [DD^] , [DDi]®p are 3x3 and [DDj] , [DD 2 ] ®p are 2x2 matrices. 

Now effective plastic stress strain matrix can be written in the 
following way. 


[D]^P 





3 j 


[DD,] 


ep 


W. 


( 2 . 20 ) 


Which represents the gauss numerical integration with 7 
gauss points of the element. [DD.|]®^ and hp are also different for 
different gauss points. 

And 


[Cjep-i = + [h-2hp) [DD^]^^]Wj (2.21), 
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Fitting [D]®p and [C,]®p in place of [D] and [C,] in the 
eq{2.6), The [K]®p matrix in place of [K] ® matrix is obtained. 

The loading and unloading parts change the characteristic 
of stress strain curve and the yield stress. This aspect is also 
included and discussed in next Chapter with other aspects in 
implementation of the computer code. 
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CHAPTER 3 


IMPLEMENTATION OF FINITE ELEMENT CODE 

This chapter explains the salient aspects of numerical 
method in developing the FEM code, including element selection. 


3.1 Element Selection 

Triangular elements are easier for discretizing complex 
geometries. In literature for different types of problems 
different type of elements are suggested as they should be free 
from defects of locking and singularity. Also the performance 
of the element under different type of loadings and boundary 
conditions should be stable. "Robust triangular element" with 
suitable nodal variables satisfy the stability conditions. 

The choice of the shape functions and variables should be 
such that the Babuska-Brezzi [Zienkiewicz and Lefebvre, 1988] 
conditions are satisfied by the system and thus, stability is 
ensured. The most vital and strictly necessary condition for the 
non- singularity of the system given by eqn(2.8) is that, 




i 1 


a 


{3.1a) 



(3.1b) 


P 


n. 


2 : 1 


n 


W 


Where ng,n^ and stand for the number of variables for 
rotation, shear force and lateral displacements respectively in 
an element. If either of the two conditions are violated, the 
locking behaviour will occur in the numerical solution. 

For the dynamic analysis if an implicit time integration 
scheme is used, it is necessary to use higher order finite 
elements and a consistent mass matrix. The higher order elements 
are effective in the representation of bending behaviour. These 
elements should be employed with consistent load vector, in that 
the midside and corner nodes are subjected to their appropriate 
load contribution in the analysis. 

In the mixed finite element formulation, triangular element 
(T6/3) has six boundary nodes which are used to determine a 
quadratic variation of 0 and w [Zienkiewicz and Lefebvre, 1988] . 
Shear forces are determined by a linear field defined by three 
internal nodes . 

When standard patch tests [Zienkiewicz and 
Lefebvre, 1988 /Zienkiewicz et al.,1986] are used with this 
element, the results indicate that, the element has incipient 
locking possibility and thus it is not a satisfactory element. 
But the examination of results indicate that the satisfaction of 
the requirements can be achieved by the addition of the 6 
variables in the interior of the element. The three internal 
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nodes created have the shape functions used of the form , 
provide a new element (T6/3 B3) which passed the same patch test 
successfully in all cases [Zienkiewicz and Lefebvre, 1988] . 

This particular element is chosen for the present analysis. 
It has 30 degree of freedom, and after the condensation of shear 
force terms the number of variables reduces to 24 . 


3.2 Time Integration scheme 

The eqn(2.8) is integrated by a numerical procedure. The 
integration is based on two ideas. 

1. Static equilibrium including inertia force is achieved at 
discrete time intervals At. 

(3.2) 

2. Variation of displacement, velocity and acceleration within 
each time step is assumed. This assumption decides the accuracy 
and stability of the solution procedure. 

The idea no. 2 decides the time integration scheme. The 
choice of a scheme depends on the finite element idealization. 
However, the finite element idealization to be chosen depends in 
turn on the actual wave velocity of the medium to be analyzed. 
It follows, therefore the selection of an appropriate finite 
element idealization of a problem and the choice of effective 
integration scheme for the response solution are closely related 
and must be considered together. 

Out of many schemes "Newmark's implicit scheme, is most 
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accurate and unconditionally stable. 


3.2a Newmark Implicit Scheme 

The Newmark' s constant acceleration scheme is written as 
follows [Bathe, 1990] , 

+ [(1-5) ^U + b^*^^U]At (3.3a) 


t + h t 


U= + ^Ut + [ (-^-a) At^ 


(3.3b) 


Where a and 6 are parameters that control integration 
accuracy and stability. Newmark originally proposed an 
unconditional stable scheme, in which case 5 = 0.5 and 
Q! = 0.25. 


Using eqn(3.2) and eqn(3.3) final form can be obtained as, 


[K+a^M] + [M] (ao^CJ + a 2 ^U + (3.4a) 


-^U) - a2U - a^U 


(3.4b) 




(3.4c) 


Where, 




The implicit scheme, as can be seen from the final 
formulation requires inversion of [K+apM] matrix. For this 
purpose Cholesky decomposition technique is computationally 
efficient. 


3.2b Time step 

After deciding on time integration scheme, it is necessary 
to use a time step suitable with finite element mesh. If the 
critical wavelength in the medium is denoted by L^, then the time 
step can be decided as, 


c n 

And the effective length of a finite element should be 

Lg = c At 


Where c is wave speed and 'n' is the number of nodes per 
critical wavelength. 

The effective wavelength and corresponding time step must 
be able to represent the complete bending wave accurately and is 
chosen appropriately depending on the kind of element 
idealization and time integration scheme used. It is reported 
that at least 8 nodes per shortest wavelength are required in 
order to produce all artifacts of wave propagation CSeron,1990] . 
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However, for higher order (parabolic and cubic) continuum 
elements the time step has to be further reduced, because the 
interior nodes are stiffer than the corner nodes [Bathe ,1990]. 
In addition for T6/3 B3 , the internal 0 variables are quartic 
functions. Also in the plate elements the flexural modes also 
affect the time step. In the actual computation appropriate time 
step is to be selected. 


3.3 Implementation Of Plastic And Dynamic Code 

At any time step the governing equation of the dynamic 
problem can be defined as eqn(2.8) . Mass matrix [M] does not get 
affected by the plasticity. So for an incremental load the 
equation can be written in the elastoplastic region. 


mnbU} + = {5f} 


(3.5a) 


Again rewriting the above equation as we did for eqn(2.17) , 
the equation acquires the following form, 

[Af]®{6CJ}.+ = {6f} + (3.5b) 


Right hand side shows the sum of the original load and the 
unbalance. The unbalance arises because of the difference in 
elastic and plastic stress strain matrices . For every time step 
this equation is balanced to get an equilibrium displacement 
satisfying the plasticity in the element. Initial unbalance is 
taken to be zero and then a trial displacement vector is 
calculated. This process is continued until both side equalize. 
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The advantage of the eqn(3.5) is that the global stiffness 
matrix is not calculated repeatedly, and thus saving computer 
time . 


3.4 Loading And Unloading 

If a material becomes plastic while loading and then 
unloaded, it behaves purely elastically. To implement this aspect 
the history of the load is to be traced to keep track of the 
strains in the material. Once the material is unloaded the yield 
point also changes. So in case of plastic loading and unloading 
the yield point of every gauss point is different. In case of 
dynamic loading the arrival of plasticity may have two effects 
on the stability. 

1. The displacements are higher than those in elastic case and 
thus functional aspect of the design gets affected. 

2. The plasticity absorbs much energy, thus the impact of 
dynamic loading may reduce. Also in case of loading and unloading 
the energy which comes under hysteresis loop gets released, so 
the damping and stability in the system may increase. 


3.5 Conqputer Implementation 

The computer code developed by Dindore [1994] has been 
modified to take the elastoplastic nature of the material in the 
present work. 

The present work deals with two type of loadings, uniform 
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distributed loading and moment loading on the boundary. Plastic 
code first determines the stresses and if these are still in the 
elastic region, the plastic part is bypassed. Also if ultimate 
stress is reached, the programme stops. If the load is in plastic 
region the load above the yield point is divided into 20 equal 
part and then applied in small steps. The core size of elasticity 
(h ) and stresses are updated after each force iteration to 
Jch computer time. '.The maximum error due to such scheme can be 

justified in following way. 

Maximum stress generated < 


so, stress increment < 


l.SOy Oy ^ 0.025O, 


20 


so the percentage of elastic height that will be missed due 
to above mentioned assumption in each force iteration is. 


Sjj < = 2.5% 


This error gets compensated to the extent because of 

following reasons, ^ vipi-iaht 

1. out Of plane stresses are not varying with heig . 

2 . weightages of various gauss points average out 

Binror" . 

3 . averaging out if the nu^er of elements are more. 

out of seven gauss points any one or more may be elastic 
Vi in the thickness. The results are 

or plastic upto any eg plasticity on the plate 

Shown as contour -plots to show extent of plastic 
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CHAPTER 4 


RESULTS AND DISCUSSIONS 


This Chapter is divided into two parts. Former consists of 
validation of computer code while later analyzes the 
interdependence of the plastic and dynamic behaviour. 


4 . 1 Validation of code 

Validation of computer code is done on following cases, 

case 1. Stresses and deflections for the purely elastic problem 
are calculated for standard cases (Timeshenko, 1959) . A 
satisfactorily agreement (within 1.5 %) is observed with 
theoretical results. Mesh is used as in fig (4. 2) . 

case 2. For qualitative validation, the plate is subjected to 
different loads and boundary conditions. The spread of 
plastic zone was observed and the same was more in the 
cases where either geometric curvature was more or there 
was a discontinuity in boundary conditions. 

case 3. For quantitative validation, a standard problem is taken 
up (Hinton, 1984) . In this a plate fixed in all edges, 
subjected to a uniform load is considered (fig 4.1) . The 
material and geometric properties are as follows, 
Modulus of elasticity, E = 30000 MPa 

Poisson's ratio, v = 0.30 

Density P = 0.003 x 10* Kg/m^ 



Thickness of plate, h = 0.20 m 

Side of the square plate, L = 6.0 m 

slope in plastic region, K = 0.01 x E 

Yield stress, 0 ^= 30.0 MPa 

Uniformly distributed load for q = 0.30, 0.35 ,0.40 MN/m^ 

The present finite element analysis consists of dividing 
the square plate by a uniform mesh as shown in fig(4.2) . As can 
be seen quarter plate is also used by applying symmetric 
conditions at the boundaries . 

Fig (4. 3) presents the deflection at the centre of the plate 
as a function of load. The figure shows the results of 
considering the mesh as a full plate and as quarter plate (which 
is equivalent to a finer mesh with elements of half size) 
together with those semi-loof analysis [Hinton and Owen, 1980] 
and Hinton and Owen, 1984. It can be seen, that the present finite 
element results are in close agreement with semi-loof analysis. 
Fig (4. 4 -4. 5) show the spread of plastic zone as contours plots 
at various load levels. At each load level the different contours 
correspond to percentage of plastic height. These plots also 
agree well with the results of Hinton (fig 4.4a) except toward 
peak load (q=0.40) Hinton's results behave more plastically. This 
happens because of layered approach used by Hinton. An inherent 
error of 12.5 percent is induced in dividing the thickness into 
8 layers. And at peak loads the present results agree with those 
of Hinton's in this error range. 

Thus, it can be seen that the results by the present code 
well with the availsble quasi-static elastoplastic 

results . 



4.2 Plastic and dynamic behaviour 


This section presents elastoplastic dynamic results 
obtained by the present approach. The variation of dynamic load 
is assumed as follows, 



COSl 



t ^ RT 


(4.1) 


For load # 1, 


Q = Qo t > RT 


And for load # 2, 


g = 


Qo cos 

\ 


Tt ( t - RT) \ 
2 RT I 


t > RT 


Where RT is the rising time for the load to reach the peak. 
The forcing functions are shown in fig{4.6). The element size 
and time steps are choosen in such a way so that the solution 
remains stable. 

To study the interdependence of plastic and dynamic 
behaviour, the analysis is divided into three parts. 

I. The amount of plasticity which is the ratio of effective 
stress to yield stress in the most plastic gauss point, is shown 
against time for various rise time (fig 4.7) . The load no. 1 is 
applied for this case. This is compared with the amount of 
plasticity in static case and also with onset of plasticity in 
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that elament . For this cas6 tha ainount of plasticity incrsases 
as the rise time is decreased while keeping other parameters 
same. By reducing the rise time the inertial forces play a 
significant role in inducing more plasticity. 

II. For a given rise time the dynamic response is analyzed for 
purely elastic and elastoplastic cases. The load no. 2 is applied 
in this case. The deflection is also seen for different loads 
with same rising time. The deflection at the centre is plotted 
against time steps. Initial maximum deflation is more than the 
purely elastic case. But at latter time the delection is less 
(fig 4.8). For this case the dynamic response in elastoplastic 
case is more damped than the purely elastic case. The response 
becomes more dampened as the plasticity increases in the plate. 
The dissipation of energy in form of hysteresis losses accounts 
for this damping. 

III. The number of gauss points which have gone plastic, are 
plotted against the time for different loads and given rise time. 
The load no. 2 is applied in this case (fig 4.9). For this case 
because of strain hardening the number of elastic gauss points 
is increasing during unloading. And at the next peak load the 
number of plastic gauss points is reducing than the previous one. 
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1 - static case 

2- Onset of plasticity 

3- • For rise time = 0.001 

4- For rise time = 0.01 


Fig (4. 7a): Effect of dynamic loading on the plasticity. 
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3- For rise time =0.001 

4 - ~ For rise time =0.01 


Fig (4. 7b): Effect of dynamic loading on plasticity 
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CHAPTER 5 


CONCLUSIONS AND SCOPE FOR FUTURE WORK 

5 . 1 Conclusions 

Based on the present finite elastoplastic and dynamic 
analysis of a plate under dynamic and static loading, following 
conclusions can be verified, 

1 . In case of dynamic loading on the plate the plasticity is 
more than it was in purely static load with the same peak value. 

2. The plasticity keeps on increasing as the rise time is 
decreased. And for a sufficient longer rise time there is no 
additional increase in the amount of plasticity over the static 
case . 

3 . As the plasticity is induced in the system, the energy is 
dissipated which acts as a dampener on the dynamic response and 
thus the oscillations reduces or vanishes. 

4. The deflection during first peak load is more in 
elastoplastic case than purely elastic case. And in latter load 
peaks because of damping the deflection comes down for 
elstoplastic case. 

5 . 2 Future Scope Of Work 

The mesh can be made finer and the time step can be made 
smaller to observe the sharper impact loading and plasticity more 

accurately. centr»l library 
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Several other loadings and different geometrical 
configurations can be tried. Also the cases of cyclic loading can 
be studied to see the permanent deflection after a given number 
of cycles. 

Normally the small cracks are always present in the 
material, but that can not be modelled in the finite element. So 
to model the cracks we can make use of statistical method with 
varying yield stresses at different gauss points. 

The effect of plasticity acts as an dampener in the dynamic 
response. The study can be useful where the system vibrations are 
to be avoided at the cost of permanent damage of some parts of 
the system. 

For dynamic loading the time step can be selected in a 
optimum manner so that the computational time is less with 
considerable accuracy and stability. 
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